Load libraries and settings
HVF analysis
hvf_res <- list()
for (i in 1:length(annot)) {
obj <- readRDS(sprintf("%s/%s_vb.rds", data_dir, names(annot[[i]])))
if (mode == "epsilon") {
hvf_res[[annot[[i]]]] <- scmet_hvf(scmet_obj = obj, delta_e = 0.9, delta_g = NULL,
evidence_thresh = 0.8, efdr = 0.1)
} else {
hvf_res[[annot[[i]]]] <- scmet_hvf(scmet_obj = obj, delta_e = NULL, delta_g = 0.5,
evidence_thresh = 0.8, efdr = 0.1)
}
}
For HVF detection task:
Evidence probability threshold chosen via EFDR valibration is too low.
Probability threshold set automatically equal to 'evidence_thresh'.
1078 Features classified as highly variable using:
- The 90 percentile of variable genes
- Evidence threshold = 0.8
- EFDR = 4.88%
- EFNR = 5.8%
For HVF detection task:
Evidence probability threshold chosen via EFDR valibration is too low.
Probability threshold set automatically equal to 'evidence_thresh'.
1811 Features classified as highly variable using:
- The 90 percentile of variable genes
- Evidence threshold = 0.8
- EFDR = 5.42%
- EFNR = 6.39%
For HVF detection task:
Evidence probability threshold chosen via EFDR valibration is too low.
Probability threshold set automatically equal to 'evidence_thresh'.
965 Features classified as highly variable using:
- The 90 percentile of variable genes
- Evidence threshold = 0.8
- EFDR = 3.08%
- EFNR = 3.48%
Summary of HVF analysis
Obtain summary of HVF analysis and also store results locally
Mean overdispersion relationship
Below we show the mean-overdispersion relationship together with the nonlinear trend fitted by scMET, from which we can extract the residual overdispersion estimates.



Mean - overdispersion plots coloured by HVF
Below we plot the mean - overdispersion relationship for each genomic context and we colour features as being highly variable (red) or other (grey).
set.seed(123)
for (an in 1:length(annot)) {
##############
## Add labels
##############
hv <- copy(hvf_res[[annot[[an]]]]$hvf$summary)
hv <- hv %>% as.data.table %>% setnames(c("feature_name"), c("id")) %>%
.[, anno := names(annot[[an]])]
# Map features to nearest genes
hv <- hv %>% merge(feat2genes[,c("id","anno","gene")], by = c("id","anno"))
# Store final labels
lab <- hv[is_variable == TRUE & gene %in% marker_genes$gene_name] %>%
.[,.SD[which.max(tail_prob)], by = "gene"] %>%
setorder(-epsilon) %>%
head(n = 15)
# Total number of hits
sig_hits <- hv[is_variable == TRUE, .N]
nonsig_hits <- hv[is_variable == FALSE, .N]
ylim <- max(hv$gamma, na.rm = TRUE) + 0.1
gg <- scmet_plot_mean_var(obj = hvf_res[[annot[[an]]]], y = "gamma",
task = "hvf", nfeatures = 2500) +
theme(legend.position = c(0.9, 0.92)) +
annotate("text", x = 0.17, y = ylim, size = 4,
label = sprintf("N=%d (%.1f%%)",
sig_hits, 100*(sig_hits/nonsig_hits))) +
geom_text_repel(force = 15, data = lab, aes_string(label = "gene"),
size = 4, col = "black", segment.color = "black",
segment.size = 0.3, segment.alpha = 0.35,
box.padding = unit(0.5,"lines"), show.legend = FALSE
)
print(gg)
pdf(sprintf("%s/ecker_hvf_%s_%s.pdf", out_dir, names(annot[[an]]),
mode), width = 6, height = 4.6)
print(gg)
dev.off()
}



Enrichment test of HVFs being marker genes
set.seed(123)
gg_list <- list()
df_enrich_list <- list()
for (i in 1:length(annot)) {
hv <- hvf_summary[[annot[[i]]]] %>%
merge(feat2genes[,c("id","anno","gene")], by = c("id","anno"))
scmet_enrich <- hv[is_variable == TRUE &
gene %in% marker_genes$gene_name] %>% .[, .N]
total_hits <- hv[is_variable == TRUE] %>% .[, .N]
null_distr <- c()
dt_copy <- copy(hv)
for (iter in 1:1000) {
idx <- sample(NROW(hv), total_hits)
null_distr[iter] <- dt_copy[, ran_hvf := FALSE] %>%
.[idx, ran_hvf := TRUE] %>%
.[ran_hvf == TRUE & gene %in% marker_genes$gene_name] %>% .[, .N]
}
df_enrich <- data.frame(x = null_distr)
df_enrich_list[[i]] <- df_enrich
gg <- ggplot(df_enrich, aes(x = x)) +
geom_histogram(aes(y = ..density..), position = "identity",
fill = "grey75", color = "grey65", bins = 20) +
geom_vline(xintercept = scmet_enrich, color = "#E69F00",
linetype = "dashed", size = 2) +
xlab("Enrichment null distribution") +
ggtitle(annot[[i]]) + ylab("Density") +
theme_classic() +
theme(
legend.position = c(0.85, 0.167),
plot.title = element_text(hjust = 0.5),
axis.text = element_text(color = "black", size = rel(0.8)),
axis.title = element_text(color = "black", size = rel(1.1))
)
gg_list[[i]] <- gg
}
print(cowplot::plot_grid(plotlist = gg_list, ncol = 3))
pdf(sprintf("%s/ecker_enrichment_hvf_%s.pdf", out_dir, mode),
width = 11.5, height = 4.5, useDingbats = FALSE)
quartz_off_screen
2

EFDR grid search plots

Mean methylation versus HVF tail probability
Below we show plots of the relationship between mean methylation and HVF tail probability, that shows the evidence of a feature being called as HVF. In general we do not observe any association, which implies that mean methylation is not a confounder for our HVF calling strategy.
set.seed(123)
p <- annot %>%
map(~ scmet_plot_vf_tail_prob(obj = hvf_res[[.]], x = "mu",
task = "hvf", nfeatures = 2000) +
ggtitle(.) + theme(legend.position = "top")) %>%
cowplot::plot_grid(plotlist = ., ncol = 3)
p
pdf(paste0(out_dir,"ecker_tailprob_vs_mu_", mode, ".pdf"),
width = 10, height = 3, useDingbats = FALSE)
quartz_off_screen
2

Summary plots
Proportion of features with \(\gamma >\) 0.5
to_plot <- rbindlist(hvf_summary) %>% .[,.(N = sum(gamma > 0.5) / .N), by = "anno_name"]
p <- ggbarplot(to_plot, x = "anno_name", y = "N", fill = "gray70") +
labs(x = "", y = expression(paste("Proportion of features with ", gamma, " > 0.5"))) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
p
pdf(paste0(out_dir, "ecker_hist_gamma.pdf"), width = 5, height = 3)
quartz_off_screen
2

Variability distribution across all genomic contexts
Residual overdispersion \(\epsilon\)
to_plot <- rbindlist(hvf_summary)
p <- ggdensity(to_plot, x = "epsilon", y = "..density..", fill = "anno_name", alpha = 0.3) +
scale_fill_brewer(palette = "Set1") +
labs(x = expression(paste("Residual overdispersion ", epsilon)), y = "Density") +
theme(legend.title = element_blank()) +
xlim(c(-2.7, 2.7))
p
pdf(sprintf("%s/ecker_density_epsilon.pdf", out_dir), width = 5,
height = 3, useDingbats = FALSE)
quartz_off_screen
2

Overdispersion \(\gamma\)
p <- ggdensity(to_plot, x = "gamma", y = "..density..", fill = "anno_name", alpha = 0.3) +
scale_fill_brewer(palette = "Set1") +
labs(x = expression(paste("Overdispersion ", gamma)), y = "Density") +
geom_vline(xintercept = 0.5, color = "black", linetype = "dashed", size = 2) +
theme(legend.title = element_blank())
p
pdf(sprintf("%s/ecker_density_gamma.pdf", out_dir), width = 5,
height = 3, useDingbats = FALSE)
quartz_off_screen
2

Mean methylation rate \(\mu\)
p <- ggdensity(to_plot, x = "mu", y = "..density..", fill = "anno_name", alpha = 0.3) +
scale_fill_brewer(palette = "Set1") +
labs(x = expression(paste("Mean methylation ", mu)), y = "Density") +
theme( legend.title = element_blank())
p
pdf(sprintf("%s/ecker_density_mu.pdf", out_dir), width = 5,
height = 3, useDingbats = FALSE)
quartz_off_screen
2

LS0tCnRpdGxlOiAiSFZGIGNhbGxpbmcgb24gcmVzaWR1YWwgb3ZlcmRpc3BlcnNpb24iCmF1dGhvcjogIkMuQS5LYXBvdXJhbmkgJiBSLiBBcmdlbGFndWV0IgpvdXRwdXQ6IAogIGh0bWxfbm90ZWJvb2s6CiAgICBkZl9wcmludDogcGFnZWQKICAgIGhpZ2hsaWdodDogaGFkZG9jawogICAgbnVtYmVyX3NlY3Rpb25zOiB5ZXMKICAgIHRoZW1lOiBjZXJ1bGVhbgogICAgdG9jOiB5ZXMKLS0tCgojIExvYWQgbGlicmFyaWVzIGFuZCBzZXR0aW5ncwpgYGB7cn0Kc3VwcHJlc3NQYWNrYWdlU3RhcnR1cE1lc3NhZ2VzKGxpYnJhcnkoc2NNRVQpKQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeShkYXRhLnRhYmxlKSkKc3VwcHJlc3NQYWNrYWdlU3RhcnR1cE1lc3NhZ2VzKGxpYnJhcnkocHVycnIpKQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMobGlicmFyeShnZ3Bsb3QyKSkKc3VwcHJlc3NQYWNrYWdlU3RhcnR1cE1lc3NhZ2VzKGxpYnJhcnkoZ2dyZXBlbCkpCnN1cHByZXNzUGFja2FnZVN0YXJ0dXBNZXNzYWdlcyhsaWJyYXJ5KGdncHVicikpCgojIyBJL08Kc291cmNlKCIuLi8uLi9sb2FkX3NldHRpbmdzLlIiKQpkYXRhX2RpciA8LSAifi9kYXRhc2V0cy9zY01FVF9tcy9lY2tlcjIwMTcvYWxsX2NlbGxzL2RhdGEvIgptYXJrZXJfZ2VuZXNfZmlsZSA8LSAifi9kYXRhc2V0cy9zY01FVF9tcy9lY2tlcjIwMTcvbWV0YWRhdGEvZWNrZXIyMDE3X21hcmtlcl9nZW5lcy5jc3YiCm91dF9kaXIgPC0gIn4vZGF0YXNldHMvc2NNRVRfbXMvZWNrZXIyMDE3L2FsbF9jZWxscy9odmYvIgpoaXRzX2RpciA8LSAifi9kYXRhc2V0cy9zY01FVF9tcy9lY2tlcjIwMTcvYWxsX2NlbGxzL2h2Zi9oaXRzLyIKaWYgKCFkaXIuZXhpc3RzKG91dF9kaXIpKSB7IGRpci5jcmVhdGUob3V0X2RpciwgcmVjdXJzaXZlID0gVFJVRSkgfQppZiAoIWRpci5leGlzdHMoaGl0c19kaXIpKSB7IGRpci5jcmVhdGUoaGl0c19kaXIsIHJlY3Vyc2l2ZSA9IEZBTFNFKSB9CgojIyBPcHRpb25zCiMgUGVyZm9ybSBIVkYgY2FsbGluZyBlaXRociB1c2luZyByZXNpZHVhbCBvdmVyZGlzcGVyc2lvbiAoZXBzaWxvbikKIyBvciBvdmVyZGlzcGVyc2lvbiAoZ2FtbWEpCm1vZGUgPC0gImVwc2lsb24iCiMgRGVmaW5lIGdlbm9taWMgY29udGV4dHMKYW5ub3QgPC0gbGlzdCgKICBjKCJkaXN0YWxfSDNLMjdhY19jb3J0ZXgiID0gIkRpc3RhbCBIM0syN2FjIiksCiAgYygiSDNLNG1lMV9jb3J0ZXgiID0gIkgzSzRtZTEiKSwKICBjKCJwcm9tXzIwMDBfMjAwMCIgPSAiUHJvbW90ZXJzIikKKQpgYGAKCiMgTG9hZCBtZXRhZGF0YSBpbmZvcm1hdGlvbgpgYGB7cn0KIyBNYXJrZXIgZ2VuZXMKbWFya2VyX2dlbmVzIDwtIGZyZWFkKG1hcmtlcl9nZW5lc19maWxlKQojIE1hcHBpbmcgZmVhdHVyZXMgdG8gbmVhcmVzdCBnZW5lcwpmZWF0MmdlbmVzIDwtIGZyZWFkKGlvJGZlYXR1cmVzMmdlbmVzKQpgYGAKCgojIEhWRiBhbmFseXNpcwpgYGB7cn0KaHZmX3JlcyA8LSBsaXN0KCkKZm9yIChpIGluIDE6bGVuZ3RoKGFubm90KSkgewogIG9iaiA8LSByZWFkUkRTKHNwcmludGYoIiVzLyVzX3ZiLnJkcyIsIGRhdGFfZGlyLCBuYW1lcyhhbm5vdFtbaV1dKSkpCiAgaWYgKG1vZGUgPT0gImVwc2lsb24iKSB7CiAgICBodmZfcmVzW1thbm5vdFtbaV1dXV0gPC0gc2NtZXRfaHZmKHNjbWV0X29iaiA9IG9iaiwgZGVsdGFfZSA9IDAuOSwgZGVsdGFfZyA9IE5VTEwsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGV2aWRlbmNlX3RocmVzaCA9IDAuOCwgZWZkciA9IDAuMSkKICB9IGVsc2UgewogICAgaHZmX3Jlc1tbYW5ub3RbW2ldXV1dIDwtIHNjbWV0X2h2ZihzY21ldF9vYmogPSBvYmosIGRlbHRhX2UgPSBOVUxMLCBkZWx0YV9nID0gMC41LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBldmlkZW5jZV90aHJlc2ggPSAwLjgsIGVmZHIgPSAwLjEpCiAgfQp9CnJtKG9iaikKYGBgCgojIFN1bW1hcnkgb2YgSFZGIGFuYWx5c2lzCk9idGFpbiBzdW1tYXJ5IG9mIEhWRiBhbmFseXNpcyBhbmQgYWxzbyBzdG9yZSByZXN1bHRzIGxvY2FsbHkKYGBge3J9Cmh2Zl9zdW1tYXJ5IDwtIGxpc3QoKQpmb3IgKGkgaW4gMTpsZW5ndGgoYW5ub3QpKSB7CiAgaHZmX3N1bW1hcnlbW2Fubm90W1tpXV1dXSA8LSBodmZfcmVzW1thbm5vdFtbaV1dXV0kaHZmJHN1bW1hcnkgJT4lIGFzLmRhdGEudGFibGUgJT4lCiAgICBzZXRuYW1lcygiZmVhdHVyZV9uYW1lIiwiaWQiKSAlPiUKICAgIC5bLCBhbm5vIDo9IG5hbWVzKGFubm90W1tpXV0pXQogIGlmIChtb2RlID09ICJlcHNpbG9uIikgewogICAgaHZmX3N1bW1hcnlbW2Fubm90W1tpXV1dXSA8LSBodmZfc3VtbWFyeVtbYW5ub3RbW2ldXV1dICU+JSBzZXRvcmRlcigtdGFpbF9wcm9iLCAtZXBzaWxvbikKICB9IGVsc2UgewogICAgaHZmX3N1bW1hcnlbW2Fubm90W1tpXV1dXSA8LSBodmZfc3VtbWFyeVtbYW5ub3RbW2ldXV1dICU+JSBzZXRvcmRlcigtdGFpbF9wcm9iLCAtZ2FtbWEpCiAgfQogIGZ3cml0ZShodmZfc3VtbWFyeVtbaV1dLCAKICAgICAgICAgZmlsZSA9IHNwcmludGYoIiVzL2h2Zl8lc18lcy50eHQuZ3oiLCBoaXRzX2RpciwgbmFtZXMoYW5ub3RbW2ldXSksIG1vZGUpKQogIGh2Zl9zdW1tYXJ5W1thbm5vdFtbaV1dXV0gJT4lIC5bLCBhbm5vX25hbWUgOj0gYW5ub3RbW2ldXSBdCn0KYGBgCgoKIyBNZWFuIG92ZXJkaXNwZXJzaW9uIHJlbGF0aW9uc2hpcCAKQmVsb3cgd2Ugc2hvdyB0aGUgbWVhbi1vdmVyZGlzcGVyc2lvbiByZWxhdGlvbnNoaXAgdG9nZXRoZXIgd2l0aCB0aGUgbm9ubGluZWFyIHRyZW5kIGZpdHRlZCBieSBzY01FVCwgZnJvbSB3aGljaCB3ZSBjYW4gZXh0cmFjdCB0aGUgcmVzaWR1YWwgb3ZlcmRpc3BlcnNpb24gZXN0aW1hdGVzLgpgYGB7ciwgZWNobz1GQUxTRSwgcmVzdWx0cz0naGlkZScsIG1lc3NhZ2U9RkFMU0UsIGZpZy53aWR0aD01LCBmaWcuaGVpZ2h0PTJ9CnNldC5zZWVkKDEyMykKZm9yIChhbiBpbiAxOmxlbmd0aChhbm5vdCkpIHsKICBnZzEgPC0gc2NtZXRfcGxvdF9tZWFuX3ZhcihvYmogPSBodmZfcmVzW1thbm5vdFtbYW5dXV1dLCB5ID0gImdhbW1hIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdGFzayA9IE5VTEwsIHNob3dfZml0ID0gVFJVRSwgbmZlYXR1cmVzID0gNDUwMCkKICBnZzIgPC0gc2NtZXRfcGxvdF9tZWFuX3ZhcihvYmogPSBodmZfcmVzW1thbm5vdFtbYW5dXV1dLCB5ID0gImVwc2lsb24iLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0YXNrID0gTlVMTCwgc2hvd19maXQgPSBUUlVFLCBuZmVhdHVyZXMgPSA0NTAwKQogIHBsb3QoY293cGxvdDo6cGxvdF9ncmlkKGdnMSwgZ2cyLCBuY29sID0gMikpCiAgcGRmKHBhc3RlMChvdXRfZGlyLCAiZWNrZXJfbWVhbnZhcl8iLCBuYW1lcyhhbm5vdFtbYW5dXSksICIucGRmIiksICAKICAgICAgd2lkdGggPSA4LjUsIGhlaWdodCA9IDMsIHVzZURpbmdiYXRzID0gRkFMU0UpCiAgcGxvdChjb3dwbG90OjpwbG90X2dyaWQoZ2cxLCBnZzIsIG5jb2wgPSAyKSkKICBkZXYub2ZmKCkKfQpgYGAKCiMgTWVhbiAtIG92ZXJkaXNwZXJzaW9uIHBsb3RzIGNvbG91cmVkIGJ5IEhWRiAKQmVsb3cgd2UgcGxvdCB0aGUgbWVhbiAtIG92ZXJkaXNwZXJzaW9uIHJlbGF0aW9uc2hpcCBmb3IgZWFjaCBnZW5vbWljIGNvbnRleHQgYW5kIHdlIGNvbG91ciBmZWF0dXJlcyBhcyBiZWluZyBoaWdobHkgdmFyaWFibGUgKHJlZCkgb3Igb3RoZXIgKGdyZXkpLgpgYGB7ciwgZmlnLndpZHRoPTMuNywgZmlnLmhlaWdodD0yLjh9CnNldC5zZWVkKDEyMykKZm9yIChhbiBpbiAxOmxlbmd0aChhbm5vdCkpIHsKICAjIyMjIyMjIyMjIyMjIwogICMjIEFkZCBsYWJlbHMKICAjIyMjIyMjIyMjIyMjIwogIGh2IDwtIGNvcHkoaHZmX3Jlc1tbYW5ub3RbW2FuXV1dXSRodmYkc3VtbWFyeSkgCiAgaHYgPC0gaHYgJT4lIGFzLmRhdGEudGFibGUgJT4lIHNldG5hbWVzKGMoImZlYXR1cmVfbmFtZSIpLCBjKCJpZCIpKSAlPiUKICAgIC5bLCBhbm5vIDo9IG5hbWVzKGFubm90W1thbl1dKV0KICAjIE1hcCBmZWF0dXJlcyB0byBuZWFyZXN0IGdlbmVzCiAgaHYgPC0gaHYgJT4lIG1lcmdlKGZlYXQyZ2VuZXNbLGMoImlkIiwiYW5ubyIsImdlbmUiKV0sIGJ5ID0gYygiaWQiLCJhbm5vIikpCiAgIyBTdG9yZSBmaW5hbCBsYWJlbHMKICBsYWIgPC0gaHZbaXNfdmFyaWFibGUgPT0gVFJVRSAmIGdlbmUgJWluJSBtYXJrZXJfZ2VuZXMkZ2VuZV9uYW1lXSAlPiUKICAgICAgLlssLlNEW3doaWNoLm1heCh0YWlsX3Byb2IpXSwgYnkgPSAiZ2VuZSJdICU+JQogICAgc2V0b3JkZXIoLWVwc2lsb24pICU+JQogICAgaGVhZChuID0gMTUpCiAgCiAgIyBUb3RhbCBudW1iZXIgb2YgaGl0cwogIHNpZ19oaXRzIDwtICBodltpc192YXJpYWJsZSA9PSBUUlVFLCAuTl0KICBub25zaWdfaGl0cyA8LSAgaHZbaXNfdmFyaWFibGUgPT0gRkFMU0UsIC5OXQogIHlsaW0gPC0gbWF4KGh2JGdhbW1hLCBuYS5ybSA9IFRSVUUpICsgMC4xCiAgZ2cgPC0gc2NtZXRfcGxvdF9tZWFuX3ZhcihvYmogPSBodmZfcmVzW1thbm5vdFtbYW5dXV1dLCB5ID0gImdhbW1hIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICB0YXNrID0gImh2ZiIsIG5mZWF0dXJlcyA9IDI1MDApICsKICAgIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9IGMoMC45LCAwLjkyKSkgKwogICAgYW5ub3RhdGUoInRleHQiLCB4ID0gMC4xNywgeSA9IHlsaW0sIHNpemUgPSA0LCAKICAgICAgICAgICAgIGxhYmVsID0gc3ByaW50ZigiTj0lZCAoJS4xZiUlKSIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNpZ19oaXRzLCAxMDAqKHNpZ19oaXRzL25vbnNpZ19oaXRzKSkpICsKICAgIGdlb21fdGV4dF9yZXBlbChmb3JjZSA9IDE1LCBkYXRhID0gbGFiLCBhZXNfc3RyaW5nKGxhYmVsID0gImdlbmUiKSwKICAgICAgICAgICAgICAgICAgICBzaXplID0gNCwgY29sID0gImJsYWNrIiwgc2VnbWVudC5jb2xvciA9ICJibGFjayIsIAogICAgICAgICAgICAgICAgICAgIHNlZ21lbnQuc2l6ZSA9IDAuMywgc2VnbWVudC5hbHBoYSA9IDAuMzUsIAogICAgICAgICAgICAgICAgICAgIGJveC5wYWRkaW5nID0gdW5pdCgwLjUsImxpbmVzIiksIHNob3cubGVnZW5kID0gRkFMU0UKICApCiAgcHJpbnQoZ2cpCgogIHBkZihzcHJpbnRmKCIlcy9lY2tlcl9odmZfJXNfJXMucGRmIiwgb3V0X2RpciwgbmFtZXMoYW5ub3RbW2FuXV0pLCAKICAgICAgICAgICAgICBtb2RlKSwgd2lkdGggPSA2LCBoZWlnaHQgPSA0LjYpCiAgcHJpbnQoZ2cpCiAgZGV2Lm9mZigpCn0KYGBgCgoKIyBFbnJpY2htZW50IHRlc3Qgb2YgSFZGcyBiZWluZyBtYXJrZXIgZ2VuZXMKYGBge3IsIGZpZy53aWR0aD04LCBmaWcuaGVpZ2h0PTN9CnNldC5zZWVkKDEyMykKZ2dfbGlzdCA8LSBsaXN0KCkKZGZfZW5yaWNoX2xpc3QgPC0gbGlzdCgpCmZvciAoaSBpbiAxOmxlbmd0aChhbm5vdCkpIHsKICBodiA8LSBodmZfc3VtbWFyeVtbYW5ub3RbW2ldXV1dICU+JSAgCiAgICBtZXJnZShmZWF0MmdlbmVzWyxjKCJpZCIsImFubm8iLCJnZW5lIildLCBieSA9IGMoImlkIiwiYW5ubyIpKQogIHNjbWV0X2VucmljaCA8LSBodltpc192YXJpYWJsZSA9PSBUUlVFICYgCiAgICAgICAgICAgICAgICAgICAgICAgZ2VuZSAlaW4lIG1hcmtlcl9nZW5lcyRnZW5lX25hbWVdICU+JSAuWywgLk5dCiAgdG90YWxfaGl0cyA8LSBodltpc192YXJpYWJsZSA9PSBUUlVFXSAlPiUgLlssIC5OXQogIG51bGxfZGlzdHIgPC0gYygpCiAgZHRfY29weSA8LSBjb3B5KGh2KQogIGZvciAoaXRlciBpbiAxOjEwMDApIHsKICAgIGlkeCA8LSBzYW1wbGUoTlJPVyhodiksIHRvdGFsX2hpdHMpCiAgICBudWxsX2Rpc3RyW2l0ZXJdIDwtIGR0X2NvcHlbLCByYW5faHZmIDo9IEZBTFNFXSAlPiUgCiAgICAgIC5baWR4LCByYW5faHZmIDo9IFRSVUVdICU+JQogICAgICAuW3Jhbl9odmYgPT0gVFJVRSAmIGdlbmUgJWluJSBtYXJrZXJfZ2VuZXMkZ2VuZV9uYW1lXSAlPiUgLlssIC5OXQogIH0KICAKICBkZl9lbnJpY2ggPC0gZGF0YS5mcmFtZSh4ID0gbnVsbF9kaXN0cikKICBkZl9lbnJpY2hfbGlzdFtbaV1dIDwtIGRmX2VucmljaAogIGdnIDwtIGdncGxvdChkZl9lbnJpY2gsIGFlcyh4ID0geCkpICsKICAgIGdlb21faGlzdG9ncmFtKGFlcyh5ID0gLi5kZW5zaXR5Li4pLCBwb3NpdGlvbiA9ICJpZGVudGl0eSIsIAogICAgICAgICAgICAgICAgICAgZmlsbCA9ICJncmV5NzUiLCBjb2xvciA9ICJncmV5NjUiLCBiaW5zID0gMjApICsKICAgIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IHNjbWV0X2VucmljaCwgY29sb3IgPSAiI0U2OUYwMCIsIAogICAgICAgICAgICAgICBsaW5ldHlwZSA9ICJkYXNoZWQiLCBzaXplID0gMikgKwogICAgeGxhYigiRW5yaWNobWVudCBudWxsIGRpc3RyaWJ1dGlvbiIpICsgCiAgICBnZ3RpdGxlKGFubm90W1tpXV0pICsgeWxhYigiRGVuc2l0eSIpICsKICAgIHRoZW1lX2NsYXNzaWMoKSArCiAgICB0aGVtZSgKICAgICAgbGVnZW5kLnBvc2l0aW9uID0gYygwLjg1LCAwLjE2NyksCiAgICAgIHBsb3QudGl0bGUgPSBlbGVtZW50X3RleHQoaGp1c3QgPSAwLjUpLAogICAgICBheGlzLnRleHQgPSBlbGVtZW50X3RleHQoY29sb3IgPSAiYmxhY2siLCBzaXplID0gcmVsKDAuOCkpLAogICAgICBheGlzLnRpdGxlID0gZWxlbWVudF90ZXh0KGNvbG9yID0gImJsYWNrIiwgc2l6ZSA9IHJlbCgxLjEpKQogICAgKQogIGdnX2xpc3RbW2ldXSA8LSBnZwp9CgpwcmludChjb3dwbG90OjpwbG90X2dyaWQocGxvdGxpc3QgPSBnZ19saXN0LCBuY29sID0gMykpCnBkZihzcHJpbnRmKCIlcy9lY2tlcl9lbnJpY2htZW50X2h2Zl8lcy5wZGYiLCBvdXRfZGlyLCBtb2RlKSwgCiAgICB3aWR0aCA9IDExLjUsIGhlaWdodCA9IDQuNSwgdXNlRGluZ2JhdHMgPSBGQUxTRSkKY293cGxvdDo6cGxvdF9ncmlkKHBsb3RsaXN0ID0gZ2dfbGlzdCwgbmNvbCA9IDMpCmRldi5vZmYoKQpgYGAKCgojIEVGRFIgZ3JpZCBzZWFyY2ggcGxvdHMKYGBge3IsIGVjaG89RkFMU0UsIHJlc3VsdHM9J2hpZGUnLCBtZXNzYWdlPUZBTFNFLCBmaWcud2lkdGg9NywgZmlnLmhlaWdodD0yLjJ9CnAgPC0gYW5ub3QgJT4lIAogIG1hcCh+IHNjbWV0X3Bsb3RfZWZkcl9lZm5yX2dyaWQob2JqID0gaHZmX3Jlc1tbLl1dLCB0YXNrID0gImh2ZiIpICsgCiAgICAgICAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIpICsgZ2d0aXRsZSguKSApICU+JQogIGNvd3Bsb3Q6OnBsb3RfZ3JpZChwbG90bGlzdCA9IC4sIG5jb2wgPSAzKQpwcmludChwKQpwZGYocGFzdGUwKG91dF9kaXIsICJlY2tlcl9lZmRyXyIsIG1vZGUsICIucGRmIiksIHdpZHRoID0gMTMsIGhlaWdodCA9IDMuNSkKcHJpbnQocCkKZGV2Lm9mZigpCmBgYAoKIyMgTWVhbiBtZXRoeWxhdGlvbiB2ZXJzdXMgSFZGIHRhaWwgcHJvYmFiaWxpdHkKQmVsb3cgd2Ugc2hvdyBwbG90cyBvZiB0aGUgcmVsYXRpb25zaGlwIGJldHdlZW4gbWVhbiBtZXRoeWxhdGlvbiBhbmQgSFZGIHRhaWwgcHJvYmFiaWxpdHksIHRoYXQgc2hvd3MgdGhlIGV2aWRlbmNlIG9mIGEgZmVhdHVyZSBiZWluZyBjYWxsZWQgYXMgSFZGLiBJbiBnZW5lcmFsIHdlIGRvIG5vdCBvYnNlcnZlIGFueSBhc3NvY2lhdGlvbiwgd2hpY2ggaW1wbGllcyB0aGF0IG1lYW4gbWV0aHlsYXRpb24gaXMgbm90IGEgY29uZm91bmRlciBmb3Igb3VyIEhWRiBjYWxsaW5nIHN0cmF0ZWd5LgpgYGB7ciwgZmlnLndpZHRoPTcsIGZpZy5oZWlnaHQ9Mi4xfQpzZXQuc2VlZCgxMjMpCnAgPC0gYW5ub3QgJT4lIAogIG1hcCh+IHNjbWV0X3Bsb3RfdmZfdGFpbF9wcm9iKG9iaiA9IGh2Zl9yZXNbWy5dXSwgeCA9ICJtdSIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRhc2sgPSAiaHZmIiwgbmZlYXR1cmVzID0gMjAwMCkgKyAKICAgICAgICBnZ3RpdGxlKC4pICsgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIpKSAlPiUKICBjb3dwbG90OjpwbG90X2dyaWQocGxvdGxpc3QgPSAuLCBuY29sID0gMykKcApwZGYocGFzdGUwKG91dF9kaXIsImVja2VyX3RhaWxwcm9iX3ZzX211XyIsIG1vZGUsICIucGRmIiksIAogICAgd2lkdGggPSAxMCwgaGVpZ2h0ID0gMywgdXNlRGluZ2JhdHMgPSBGQUxTRSkKcHJpbnQocCkKZGV2Lm9mZigpCmBgYAoKIyBTdW1tYXJ5IHBsb3RzCgojIyBQcm9wb3J0aW9uIG9mIGZlYXR1cmVzIHdpdGggJFxnYW1tYSA+JCAwLjUKYGBge3J9CnRvX3Bsb3QgPC0gcmJpbmRsaXN0KGh2Zl9zdW1tYXJ5KSAlPiUgLlssLihOID0gc3VtKGdhbW1hID4gMC41KSAvIC5OKSwgYnkgPSAiYW5ub19uYW1lIl0KcCA8LSBnZ2JhcnBsb3QodG9fcGxvdCwgeCA9ICJhbm5vX25hbWUiLCB5ID0gIk4iLCBmaWxsID0gImdyYXk3MCIpICsKICBsYWJzKHggPSAiIiwgeSA9IGV4cHJlc3Npb24ocGFzdGUoIlByb3BvcnRpb24gb2YgZmVhdHVyZXMgd2l0aCAiLCBnYW1tYSwgIiA+IDAuNSIpKSkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gMCwgaGp1c3QgPSAwLjUpKQpwCnBkZihwYXN0ZTAob3V0X2RpciwgImVja2VyX2hpc3RfZ2FtbWEucGRmIiksIHdpZHRoID0gNSwgaGVpZ2h0ID0gMykKcHJpbnQocCkKZGV2Lm9mZigpCmBgYAoKIyMgVmFyaWFiaWxpdHkgZGlzdHJpYnV0aW9uIGFjcm9zcyBhbGwgZ2Vub21pYyBjb250ZXh0cwoKIyMjIFJlc2lkdWFsIG92ZXJkaXNwZXJzaW9uICRcZXBzaWxvbiQKYGBge3J9CnRvX3Bsb3QgPC0gcmJpbmRsaXN0KGh2Zl9zdW1tYXJ5KQpwIDwtIGdnZGVuc2l0eSh0b19wbG90LCB4ID0gImVwc2lsb24iLCB5ID0gIi4uZGVuc2l0eS4uIiwgZmlsbCA9ICJhbm5vX25hbWUiLCBhbHBoYSA9IDAuMykgKwogIHNjYWxlX2ZpbGxfYnJld2VyKHBhbGV0dGUgPSAiU2V0MSIpICsKICBsYWJzKHggPSBleHByZXNzaW9uKHBhc3RlKCJSZXNpZHVhbCBvdmVyZGlzcGVyc2lvbiAiLCBlcHNpbG9uKSksIHkgPSAiRGVuc2l0eSIpICsKICB0aGVtZShsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpICsKICB4bGltKGMoLTIuNywgMi43KSkKcApwZGYoc3ByaW50ZigiJXMvZWNrZXJfZGVuc2l0eV9lcHNpbG9uLnBkZiIsIG91dF9kaXIpLCB3aWR0aCA9IDUsIAogICAgaGVpZ2h0ID0gMywgdXNlRGluZ2JhdHMgPSBGQUxTRSkKcHJpbnQocCkKZGV2Lm9mZigpCmBgYAoKIyMjIE92ZXJkaXNwZXJzaW9uICRcZ2FtbWEkCmBgYHtyfQpwIDwtIGdnZGVuc2l0eSh0b19wbG90LCB4ID0gImdhbW1hIiwgeSA9ICIuLmRlbnNpdHkuLiIsIGZpbGwgPSAiYW5ub19uYW1lIiwgYWxwaGEgPSAwLjMpICsKICBzY2FsZV9maWxsX2JyZXdlcihwYWxldHRlID0gIlNldDEiKSArCiAgbGFicyh4ID0gZXhwcmVzc2lvbihwYXN0ZSgiT3ZlcmRpc3BlcnNpb24gIiwgZ2FtbWEpKSwgeSA9ICJEZW5zaXR5IikgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAuNSwgY29sb3IgPSAiYmxhY2siLCBsaW5ldHlwZSA9ICJkYXNoZWQiLCBzaXplID0gMikgKwogIHRoZW1lKGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKcApwZGYoc3ByaW50ZigiJXMvZWNrZXJfZGVuc2l0eV9nYW1tYS5wZGYiLCBvdXRfZGlyKSwgd2lkdGggPSA1LCAKICAgIGhlaWdodCA9IDMsIHVzZURpbmdiYXRzID0gRkFMU0UpCnByaW50KHApCmRldi5vZmYoKQpgYGAKCgojIyMgTWVhbiBtZXRoeWxhdGlvbiByYXRlICRcbXUkCmBgYHtyfQpwIDwtIGdnZGVuc2l0eSh0b19wbG90LCB4ID0gIm11IiwgeSA9ICIuLmRlbnNpdHkuLiIsIGZpbGwgPSAiYW5ub19uYW1lIiwgYWxwaGEgPSAwLjMpICsKICBzY2FsZV9maWxsX2JyZXdlcihwYWxldHRlID0gIlNldDEiKSArCiAgbGFicyh4ID0gZXhwcmVzc2lvbihwYXN0ZSgiTWVhbiBtZXRoeWxhdGlvbiAiLCBtdSkpLCB5ID0gIkRlbnNpdHkiKSArCiAgdGhlbWUoIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKcApwZGYoc3ByaW50ZigiJXMvZWNrZXJfZGVuc2l0eV9tdS5wZGYiLCBvdXRfZGlyKSwgd2lkdGggPSA1LCAKICAgIGhlaWdodCA9IDMsIHVzZURpbmdiYXRzID0gRkFMU0UpCnByaW50KHApCmRldi5vZmYoKQpgYGAKCg==